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In order to obtain an evolution system which is robust against the violation of constraints, we 
present a new set of evolution systems based on the so-called Baumgarte-Shapiro-Shibata-Nakamura 
(BSSN) equations. The idea is to add functional derivatives of the norm of constraints, C 2 , to the 
evolution equations, which was proposed by Fiske (2004) and was applied to the ADM formulation 
in our previous study. We derive the constraint propagation equations, discuss the behavior of 
constraint damping, and present the results of numerical tests using the gauge-wave and polarized 
Gowdy wave spacetimes. The construction of the C 2 -adjusted system is straightforward. However, 
in BSSN, there are two kinetic constraints and three algebraic constraints; thus, the definition of 
C 2 is a matter of concern. By analyzing constraint propagation equations, we conclude that C 2 
should include all the constraints, which is also confirmed numerically. By tuning the parameters, 
the lifetime of the simulations can be increased as 2-10 times as longer than those of the standard 
BSSN evolutions. 

PACS numbers: 04.25.D- 



I. INTRODUCTION 

When solving the Einstein equations numerically, the 
standard way is to split the spacetime into space and 
time. The most fundamental decomposition of the Ein- 
stein equations is the Arnowitt-Deser-Misner (ADM) for- 
mulation 0, 0). However, it is well known that in long- 
term evolutions in strong gravitational fields such as the 
coalescences of binary neutron stars and/or black holes, 
simulations with the ADM formulation are unstable and 
are often interrupted before producing physically inter- 
esting results. Finding more robust and stable formula- 
tions is known to the "formulation problem" in numerical 
relativity 0-11]. 

Many formulations have been proposed in the last two 
decades. The most commonly used sets of evolution 
equations among numerical relativists are the so-called 
Baunrgarte-Shapiro-Shibata-Nakamura (BSSN) formula- 
tion Q 0], the generalized harmonic (GH) formulation 
, the Kidder-Scheel-Teukolsk y ( KST) formulation 
and the Z4 formulation [ill . Tl2| (as references of 
their numerical application, we here cite only well-known 
articles; [H, Q3] for the BSSN formulation, [H[ for the 
GH formulation, |l6| for the KST formulation, and 17 1 
for the Z4 formulation). 

All of the above modern formulations include the tech- 
nique of "constraint damping" , which attempts to con- 
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trol the violations of constraints by adding the constraint 
terms to their evolution equations. Using this technique, 
more stable and accurate systems are obtained (see e.g. 
OH!!])- This technique can be described as 'adjustment' 
of the original system. 

In [20H22j . two of the authors systematically investi- 
gated how the adjusted terms change the original sys- 
tems by calculating the constraint propagation equations. 
The authors suggested some effective adjustments for 
the BSSN formulation under the name "adjusted BSSN 
formulation" [22j . The actual constraint-damping effect 
was confirmed by numerical tests (23j . 

Fiske proposed a method of adjusting the original evo- 
lution system using the norm of the constraints, C 2 , [24| . 
which wc call a "C 2 -adjusted system." The new evo- 
lution equations force the constraints to evolve towards 
their decay if the coefficient parameters of the adjusted 
terms are set as appropriate positive values. Fiske re- 
ported the damping effect of the constraint violations for 
the Maxwell system I24f and for the linearized ADM and 
BSSN formulations |25|. He also reported the limitation 
of the magnitude of the coefficient parameters of the ad- 
justed terms. 

In (26|, we applied this C 2 -adjusted system to the (full) 
ADM formulation and presented some numerical tests. 
We confirmed that the violations of the constraints are 
less than those in the original system. We also reported 
the differences of the effective range of the coefficient of 
the adjusted terms. 

In this article, we apply the C 2 -adjusted system to the 
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(full) BSSN formulation and derive the constraint propa- 
gation equations in the flat space. We perform some nu- 
merical tests and compare them with three other types 
of BSSN formulations: the standard BSSN formulation, 
the A-adjusted BSSN formulation, and the C 2 -adjustcd 
BSSN formulation. We use the gauge- wave and polarized 
Gowdy wave testbeds, which are the test problems as is 
known to apples- with-apples testbeds for comparing evo- 
lution systems [I?}. Since the models are precisely fixed 
up to the gauge conditions, boundary conditions, and 
technical parameters, the testbeds are widely used for 
comparisons [H H, [29| . 

The structure of this article is as follows. We review 
the ideas of adjusted systems and C 2 -adjusted system in 
Scc|n] In Scc lIII| we review the standard and adjusted 
BSSN formulations and derive the C 2 -adjusted version 
of the BSSN formulation. In Scc lIV| we present some 
numerical tests of the gauge-wave and polarized Gowdy 
wave testbeds. We show the damping effect of the con- 
straint violations, and confirm that inclusion of algebraic 
constraints in C 2 make the violations of constraints de- 
crease. We summarize this article in Sec|V] In this arti- 
cle, we only consider vacuum spacetime, but the inclusion 
of matter is straightforward. 

II. IDEAS OF ADJUSTED SYSTEMS AND 
C* 2 -ADJUSTED SYSTEMS 

A. Idea of adjusted systems 

Suppose we have dynamical variables u % that evolve 
with the evolution equations 

8tu* ; = /(«*,«,•«',••■)> (2.1) 

and suppose also that the system has the (first class) 
constraint equations 

C'V^uY--) »0. (2.2) 

We can then predict how the constraints are preserved 
by evaluating the constraint propagation equations 

d t C a =g(C a ,d l C a , ■••), (2.3) 

which measure the violation behavior of constraints C a 
in time evolution. Equation (|2.3p is theoretically weakly 
zero, i.e., d t C a ~ 0, since the system is supposed to 
be the first class. However, free numerical evolution 
with discretized grids introduces a constraint violation, 
at least at the level of truncation error, which sometimes 
grows and stops the simulations. The unstable feature of 
ADM evolution can be understood on the basis of this 
analysis [l5| . 

Such features of the constraint propagation equations, 
(j2.3[) , change when we modify the original evolution equa- 
tions. Suppose we add constraint terms to the right- 
hand-side of (|2.1[) as 

dtu* = f(u i ,a j u i ,'--)+F{C a ,d j C a ,-"), (2.4) 



where F(C a , ■■•) « in principle zero but not exactly 
zero in numerical evolutions. With this adjustment, 
equation (|2.3p will also be modified to 

d t C a =g{C a ,d i C a ,---)+G{C a ,d i C a ,---). (2.5) 

Therefore, we are able to control dtC a by making an 
appropriate adjustment F(C°, djC a , ■ ■ ■ ) in (g^j) . If 
dtC a < is realized, then the system has the constraint 
surface as an attractor. 

This technique is also known as a constraint-damping 
technique. Almost all the current popular formulations 
used in large-scale numerical simulations include this im- 
plementation. The purpose of this article is to find a 
better way of adjusting the evolution equations to realize 
d t C a < 0. 



B. Idea of C 2 -adjusted systems 

Fiske [24| proposed a way of adjusting the evolution 
equations which we call "C 2 -adjusted systems" ; 

= /(«*, a,-**,--. , (2.6) 

where n 13 is a positive-definite constant coefficient and 
C 2 is the norm of the constraints, which is defined as 

C 2 = J C a C a d 3 x. The term (SC 2 /Su j ) is the functional 

derivative of C 2 with respect to u 3 . The associated con- 
straint propagation equation becomes 

<W-mc-,*o-,-)-/-'.(«£) 

(2.7) 

The motivation for this adjustment is to naturally ob- 
tain the constraint-damping system, dtC 2 < 0. If we set 
k 13 so that the second term of the right-hand side of (|2.7[) 
becomes larger than the first term, then dtC 2 becomes 
negative, which indicates that constraint violations are 
expected to decay to zero. Fiske presented numerical ex- 
amples of the Maxwell system and the linearized ADM 
and BSSN formulations, and concluded that this method 
actually reduces constraint violations as expected. In our 
previous work (26|, we applied the C 2 -adjusted system to 
the (full) ADM formulation and derived the constraint 
propagation equations. We confirmed that dtC 2 < is 
expected in the flat spacetime. We performed numerical 
tests with the C 2 -adjusted ADM formulation using the 
Gowdy wave testbed, and confirmed that the violations 
of the constraint are lower than those of the standard 
ADM formulation. The simulation continues 1.7 times 
longer than that of the standard ADM formulation with 
the magnitude of the violations of the constraint less than 
order 0(10°). 
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III. APPLICATION TO BSSN FORMULATION 
A. Standard BSSN Formulation 

We work with the widely used notation of the 
BSSN system. That is, the dynamical variables 
(ip, K, jij, Aij, r l ) as the replacement of the variables of 
the ADM formulation, (jij, K^j), where 



(3.1) 
(3.2) 
(3-3) 



p=(l/12)log(det(7«)). 
K -'''A,,. 

— e ^ , 

= e~ Av {K l3 - (l/3) 7l3 A'), and (3.4) 

The BSSN evolution equations are, then, 

d t tp = -(l/6)aK + (l/6)(a 4 /3 l ) + (3 l (d lV ), (3.6) 
d t K = aAijAV + (1/3)q/v 2 - D t D l a + (3 l (d l K), 

(3.7) 

dfiij = -2aAij - (2/3)7 y -(ft^) 

+ i it {dip) + luidjP*) + PidfYij), (3-8) 
d t A zj = aKAij - 2aA u A i j + ae~ A<p Rij T¥ 

- e-^(D t D ]a ) TV - {2/2,)A ij (d^) 

+ (diftAjt + (d^ e )Au + PldtAit), (3.9) 
dtf 1 = 2a{6(d j <p)A ij +fV^ - (2/3)i ij (dj K)} 

- 2(0 i a)A« + (2/3)f i (5 i ^) + (1/3)^(8^^) 

+ f3 e (d e r) - r^djp^ + ^idjdeF), (3.10) 

where TF denotes the trace-free part. The Ricci tensor 
in the BSSN system is normally calculated as 



Rij — Rij + ijy , 



(3.11) 



where 



ln(id 3 )T n + 7 £m (2r fe f ( i r j ) fcm + T n ijT n im ) 



x (iUj)L + 7 (21 j) A 

(i/2)7 m ^, m£ + f"r fo)r 



(3.12) 



Rf. = -2A-D^ + 4(A<^)(-Diy) - 2y i:i D m D m tp 



4% (£>•>)(£>„,¥>). 



(3.13) 



The BSSN system has five constraint equations. The 
"kinematic" constraint equations, which are the Hamilto- 
nian constraint equation and the momentum constraint 
equations (H-constraint and .M-constraint, hereafter), 
are expressed in terms of the BSSN basic variables as 

H = e-^R - 8e- 4v (A;5V + (D m <p)(D m ip)) 

+ (2/3) A" 2 - A tj A l i - (2/3) AK « 0, (3.14) 
Mi = -{2/3)D l K + 6(5^)1^ + DjA'i 

- 2{b iV ,)A « 0, (3.15) 



respectively, where Di is the covariant derivative associ- 
ated with 7y and R — 7 13 R^. Because of the introduc- 
tion of new variables, there are additional "algebraic" 
constraint equations: 



gi = p _ jj^j 



0. 



.4 = A«7y « 0, 

5 = det(7y) - 1 « 0, 



(3.16) 

(3.17) 
(3.18) 



which we call the A-, and 5-constraints, respectively, 
hereafter. If the algebraic constraint equations, (|3.16D - 
p,18|) . are not satisfied, the BSSN formulation and ADM 
formulation are not equivalent mathematically. 



B. C 2 -adjusted BSSN Formulation 

The C 2 -adjusted BSSN evolution equations are for 
mally written as 

w ( ,i 9 , 

(3.20) 
(3.21) 
(3.22) 
(3.23) 



1 Stp J ' 
SC 2 ' 



d t <p = msn - a 

d t K = (I3T71) - \k 

dtjij = (|3 ■ 8^ — X^ijmn 



5K J ' 
SC 2 



SC 2 



5A r 



r \6T3 

where all the coefficients X v , Xk , X^ij mn , X^ijmn > and A~ 
are positive definite. C 2 is a function of the constraints 
"H, Mi, G\ A, and 5, which we set as 



C 2 



If 2 ■ - •■M,M, ■ c <rij Q i Q j 
f c A A 2 + c s S 2 )d 3 x, 



(3.24) 

where, cq, c-a-i an d eg are Boolean parameters (0 or 1). 
These three parameters are introduced to prove the ne- 
cessity of the algebraic constraint terms in (|3.24[) . 

The adjusted terms in (|3.19[) - (|3.23p are then written 
down explicitly, as shown in Appendix |XJ The constraint 
propagation equations of this system are also derived for 
the Minkowskii background, as shown in Appendix [Bl 

Now we discuss the effect of the algebraic constraints. 
From (|B1|) - (|B5|) . we see that the constraints affect each 
others. The constraint propagation equations of the 
algebraic constraints, (|B3[) - (|B5|) . include cc(A^A<5 a t — 
2Xf5 a b)Q b ', — 6caA^^4, and — 6csXjS, respectively. These 
terms contribute to reduce the violations of each con- 
straint if cq, CA-, and cs are non-zero. Therefore, we 
adopt cq = ca — cs — 1 in (|3.24l) : 

C 2 = J (U 2 + ^MiMj + 7u W + A 2 + S 2 ) d 3 x. 

(3.25) 
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This discussion is considered only from the viewpoint of 
the inclusion of the diffusion terms. In order to validate 
this decision, we perform some numerical examples in 
SecED 



C. ,4-adjusted BSSN System 

In [22|], two of the authors reported some examples 
of adjusted systems for the BSSN formulation. The au- 
thors investigated the signatures of eigenvalues of the co- 
efficient matrix of the constraint propagation equations, 
and concluded three of the examples to be the best can- 
didates for the adjustment. The actual numerical tests 
were performed later [23[ using the gauge-wave, linear- 
wave, and polarized Gowdy wave testbeds. The most 
robust system among the three examples for these three 
testbeds was the A-adjusted BSSN formulation, which 
replaces (|3.9p in the standard BSSN system with 



KAOtD(iMj), 



(3.26) 



where ka is a constant. If ka is set as positive, the 
violations of the constraints are expected to be damped 
in flat spacetime 22j. We also use the ^-adjusted BSSN 



system for comparison in the following numerical tests. 



IV. NUMERICAL EXAMPLES 

We test the three systems (C 2 -adjusted BSSN, A- 
adjustcd BSSN, and standard BSSN) in numerical evo- 
lutions using the gauge-wave and polarized Gowdy wave 
spacetimes, which are the standard tests for comparisons 
of formulations in numerical relativity, and are known 
as apples- with-applcs testbeds [27|- We also performed 
the linear-wave testbed but the violations of the con- 
straint are negligible; thus, we employ only the above 
two testbeds in this article. These tests have been used 
by severa l group s and were reported in the same manner 

(e.g., iniiig). 

For simplicity, we set the coefficient parameters in 

(|3.2I[)-([3,23j) to Xjijmn = XjSimSjn, ^Aijmn = ^A^ im ^3 n > 



Ap<5 y with non- negative coefficient constant 



and A~ 
r 

parameters A^, Xj, and Ap. Our code passes the conver- 
gence test with second-order accuracy. We list the figures 
in this article in Table U for reader's convenience. 



A. Gauge-wave Testbed 

1. Metric and Parameters 

The metric of the gauge- wave test is 

ds 2 = -Hdt 2 + Hdx 2 + dy 2 + dz 2 , (4.1) 
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(A-1) H-constraint 
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FIG. 1. L2 norm of each constraint violation in the gauge- 
wave evolution using the standard BSSN formulation. The 
vertical axis is the logarithm of the L2 norm of the constraints 
and the horizontal axis is time. We see the evolution stops at 
t = 110 due to the growth of .M-constraint violation. 



where 



H= 1 - A sm(2ir(x -t)/d), 



(4.2) 



which describes a sinusoidal gauge wave of amplitude A 
propagating along the x-axis. The nontrivial extrinsic 
curvature is 



K x 



ttA cos( — y ~2 — -) 

I - Asm — y ~2 — - 



(4.3) 



Following [27] . we chose the numerical domain and pa- 
rameters as follows: 

• Gauge- wave parameters: d = 1 and A = I0~ 2 . 

• Simulation domain: x G [—0.5, 0.5], y — z = 0. 

• Grid: x n = -0.5 + (n-l/2)dx with n = 1, • • • ,100, 
where dx = I/I00. 

• Time step: dt = 0.25dx. 

• Boundary conditions: Periodic boundary condition 
in x-direction and planar symmetry in y- and z- 
directions. 

• Gauge conditions: 

d t a = -a 2 K, ft = 0. (4.4) 

• Scheme: second-order iterative Crank-Nicolson. 

2. Constraint Violations and Their Dampings 

Figure [T] shows the violations of five constraint equa- 
tions H, Aii, Q l , A, and S for the gauge-wave evolution 
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TABLE I. List of figures. 

gauge-wave test Gowdy wave test 

gvg aisgf 

(A) standard BSSN (I3.6p - (|3.f Op FigJTJnorm each Fig® norm each 
(constraint propagation, see App. [Cj Fig[2]norm all Fig[7]norm all 

(B) ^-adjusted BSSN Fig[2]norm all Fig[7]norm all 
{32J-J33]), (pAOl) , and ([3~26j) Fig[3]norm each 

(constraint propagation, see App. iBl) 

(C) C^-adiusted BSSN f|3. 19|) - (|3.23| i Fig[2]norm all Fig[7]norm all 
(constraint propagation, see App. [B} Fig[3]norm each Fig[8]norm each 

FigHl adjusted ratio Fig[9] adjusted ratio 
Fig[5] ^25\i test Fig[T0] (p^5|) test 
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(B) A-adjusted BSSN : 

(C) C 2 -adjusted BSSN ! 
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FIG. 2. L2 norm of all the constraints in gauge-wave 
evolution comparing three BSSN formulations: (A) standard 
BSSN formulation (solid line), (B) A-adjusted BSSN formu- 
lation (dotted line), and (C) C 2 -adjusted BSSN formulation 
(dot-dashed line). The adopted parameters are ka = 1CP 1 ' 6 
for (B), and \ v = 10" 8 ' 5 , \ K = 10" 8 ' 4 , A^ = 10" 7 ' 3 , 
A^ = 10~ 2 ' 5 , and A f = 10 -1 ' 8 for (C) to minimize C 2 at 
t = 1000. The constraint violations of the A-adjusted BSSN 
formulation, (B), increase with time and the simulation stops 
before t = 1300, while those of the C 2 -adjusted BSSN for- 
mulation, (C), remain at O(10 _1 ) until t = 1300 and the 
simulation stops at t = 1350. 



using the standard BSSN formulation. The violation of 
the A^-constraint, line (A-2), is the largest during the 
evolution, while the violations of both the .A-constraint 
and iS-constraint are negligible. This is the starting point 
for improving the BSSN formulation. 

Applying the adjustment procedure, the lifetime of the 
standard BSSN evolution is increased at least 10-fold. In 
Figl^l we plot the L2 norm of the constraints, (|3.25[) . of 
three BSSN evolutions: (A) the standard BSSN formula- 
tion (|3~B1) - (PUJ1) (B) th e A-a djusted BSSN formulation 
dSUl-jSU), dSHUJ), and ( E T261), and (C) the C 2 -adjusted 
BSSN formulation ((3~T3)) - ((3~23)) . For the standard BSSN 
case, we see the violation of constraint monotonically in- 
creases in the earlier stage, while other two adjusted cases 



keep it smaller. We can say that the C 2 -adjusted for- 
mulation is the most robust one against the violation of 
constraints between three. 

We plot the norm of each constraint equation in Fig|3l 
First, we see that the violation of the A'l-constraint for 
the two adjusted BSSN formulations [the lines (B-2) and 
(C-2) in Fig[3] are less than that of the standard BSSN 
formulation in FigJTJ This behavior would be explained 
from the constraint propagation equations, where we see 
the terms X^AM a and (l/2)n A AMi in jB2| and (fB7| . 
respectively. These terms contribute to reduce the viola- 
tions of the Al-constraint. This is the main consequence 
of the two adjusted BSSN formulations. 

Second, we also find that the violations of the A- 
constraint and 5-constraint are larger than those in FigJT] 
From constraint propagation equations (|B4[) and (|C4[) . 
the violation of the „4-constraint is triggered by the Al- 
and ,4-constraints. The increase in the violations of the 
^4-constraint is caused by the term 2A^(5 y (d-iAij). Simi- 
larly, in (|B5[) and (|C5[) , the violation of the 5-constraint 
is triggered by only the ,4-constraint since the magni- 
tude of A=y is negligible. Therefore, the increase in the 
violation of the iS-constraint is due to the violation of 
the ^-constraint. 

From (|A1[) and (|A3[) . it can be seen that the ad- 
justed terms of the evolution equations of ip and jij in- 
clude second-order derivative terms of the ^-constraint. 
This means that these evolution equations include fourth- 
order derivative terms of the dynamical variables. In or- 
der to investigate the magnitudes of the adjusted terms, 
we show in Fig[4]the ratio of the adjusted terms to that 
of the original terms in each evolution equation. We see 
that the magnitudes of the adjusted terms of ip and 7^ 
are reasonably small. 

In the simulations with the C 2 -adjusted BSSN formu- 
lation, the largest violation is the iS-constraint. The iS- 
constraint depends only on the dynamical variables 7^, 
so that there is no other choice than setting A^ for con- 
trolling 5-constraint, as can be seen from (|B5[) . How- 
ever, we must set Xj to a value as small as possible since 
the adjusted term of 7^ includes higher derivatives of 
7y. Therefore, it is hard to control the 5-constraint, and 
we have not yet found an appropriate set of parameters. 
This will remain as a future problem of this C 2 -adjusted 
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FIG. 3. L2 norm of each constraint in the gauge-wave evolution using the A-adjusted BSSN formulation [panel (a)] and 
C 2 -adjusted BSSN formulation [panel (b)]. The parameters ka, X v , Xk, X~ f , At, and Ap are the same as those in Fig[2] In 
both panels, we see that the violations of the "H-constraint [the lines (B-1) and (C-1)], the .M-constraint [(B-2) and (C-2)], and 
the ^-constraint [(B-3) and (C-3)] are less than those for the standard BSSN formulation in Fig[TJ However, the violations of 
the ^-constraint [(B-4) and (C-4)] and the iS-constraint [(B-5) and (C-5)] are larger. Line (B-5) overlaps with line (B) in Fig[5] 
after t = 100, and line (C-5) overlaps with line (C) in Fig[2] after t = 500. 
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FIG. 4. L2 norm of the ratio (adjusted terms)/(original 
terms) of each evolution equation of the C 2 -adjusted BSSN 
formulation, (|3.19[) - (|3.23[) . in the gauge-wave test. We see 
that the largest ratio is the evolution equation of Aij. The 
corrections to (p, K, and jij evolution equations are reason- 
ably small. 



BSSN system. 

We also investigated the sensitivity of the parame- 
ters in the C 2 -adjusted BSSN evolutions. We com- 
pared evolutions with setting only one of the param- 
eters, (A v , Xk, A;y, At, Ap), nonzero. Since the key of 
the damping of the violation of constraints is the A4- 
constraint, and (Air, A t) c ontrols the violation of M- 
constraint directly by (|B2[) . we mention here only the 
dependence on Xk and A^. We found that constraint- 
damping feature changes sensitively by both Xk and A^j, 
among them setting Aj^ is important to control the M.- 
constraint violation. Wc sec the best controlled evolution 



FIG. 5. Difference with the definition of C 2 , (J3T25J , in the 
damping of each constraint violation with cg = ca = cs = 0. 
The parameters \ v , Xk, Xf, At, and Ap are the same as 
those in Fig(2] The simulation stops since the violations of 
the constraints sudden increase at t = 800. 

with X A = 1CT 3 , than 1CT 2 and 1CT 4 . 



3. Contribution of Algebraic Constraints 
in Definition of C 2 

In Scc lIIIBI we defined C 2 , ([3T2g]> . including the alge- 
braic constraints. We check this validity by turning off 
the algebraic constraints in (|3.25p . The result is shown in 
FigfSJ where we see the simulation stops at t = 800 due to 
a sudden increase in the violation of the constraints. This 
confirms that the algebraic constraints play an important 
role of damping of the violations of constraints. We also 
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tested with other combinations of Boolean parameters 
(cgjCAiCs): and confirmed that the best controlled evo- 
lution is realized when cq — ca = eg = 1. 



B. Gowdy-wave Testbed 

1. Metric and Parameters 

The metric of the polarized Gowdy wave is given by 

ds 2 = r^e^i-dt 2 + dx 2 ) + t(e p dy 2 + e^dz 2 ), 

(4.5) 

where P and A are functions of x and t. The forward 
direction of the time coordinate t corresponds to the ex- 
panding universe, and t = corresponds to the cosmo- 
logical singularity. 

For simple forms of the solutions, P and A are given 
by 

P = J (2irt)cos(2TTx), (4.6) 
A = -27rfJ (27rf) Ji(27rf) cos 2 (2ttx) + 27r 2 t 2 [J 2 (2vrt) 
+ J 2 (2nt)} (l/2){(2n) 2 [J 2 (2n) + J 2 (2^)] 
-27rJ (27r)Ji(27r)}, (4.7) 

where J„ is the Bessel function. 

Following [27]], a new time coordinate r, which satisfies 
harmonic slicing, is obtained by the coordinate transfor- 
mation 

t(r) = ke CT , (4.8) 

where k and c are arbitrary constants. We also follow 
{27} by setting k, c, and the initial time to as 

k ~ 9.67076981276405, c - 0.002119511921460, 

(4.9) 

t = 9.87532058290982, (4.10) 

so that the lapse function in the new time coordinate is 
unity and t = r at the initial time. 

We also use the following parameters specified in [27| . 

• Simulation domain: x £ [—0.5, 0.5], y = z = 0. 

• Grid: x„ = -0.5 + (n - (l/2))da;, n = 1, • ■ • , 100, 
where dx = 1/100. 

• Time step: dt = 0.25dx. 

• Boundary conditions: Periodic boundary condition 
in x-direction and planar symmetry in y- and z- 
directions. 

• Gauge conditions: dta = —a 2 K, fi l = 0. 

• Scheme: second-order iterative Crank-Nicolson. 




200 400 600 800 1000 

-Time 

FIG. 6. L2 norm of each constraint equation in the polarized 
Gowdy wave evolution using the standard BSSN formulation. 
The vertical axis is the logarithm of the L2 norm of the con- 
straint and the horizontal axis is backward time. 

2. Constraint Violations and Their Dampings 

We begin showing the case of the standard BSSN for- 
mulation, (|3.6[) - (|3.10|) . Figure |B] shows the L2 norm of 
the violations of the constraints as a function of back- 
ward time (—t). We see that the violation of the A4- 
constraint is the largest at all times and that all the vi- 
olations of constraints increase monotonically with time. 
[Comparing with the result in [23|, our code shows that 
the H-constraint (A-l) remains at the same level but the 
A4-constraint (A-2) is smaller.] 

Similar to the gauge- wave test, we compare the viola- 
tions of C 2 for three types of BSSNs in FigJT] In the case 
of the A-adjusted BSSN formulation, the violation of the 
constraints increases if we set \ka\ larger than 10~ a2 . 
In the case of the C 2 -adjusted BSSN formulation, it in- 
creases if we set |Ar| larger than 10 -12 . Note that the 
signatures of the above ka and As are negative, contrary 
to the predictions in [22j and Sec lIIH respectively. This is 
because these simulations are performed with backward 
time. 

As shown in FigJTl the violations of C 2 for the standard 
BSSN formulation and the A-adjusted BSSN formulation 
increase monotonically with time, while that for the C 2 - 
adjusted BSSN formulation decreases after t = —200. 
To investigate the reason of this rapid decay after t = 
—200, we plot each constraint violation in FigJS] We see 
that the violations of the .A-constraint and 5-constraint 
increase with negative time, in contrast to the standard 
BSSN formulation, and those of the A4-constraint and 
^-constraint decrease after t = —200. The propagation 
equation of the A4-constraint, (|B2[) , includes the term 
— 2cAX^d a A, which contributes to constraint damping. 
Similarly, the propagation equation of the ^-constraint, 
(|B3jl . includes 5 ab {(l/2)\^d b A + 2\ f d b }H~c s \^5 ab d b S- 1 
the decay of the violations of the ^-constraint is caused 
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200 400 600 800 1000 200 400 600 800 1000 



-Time -Time 



FIG. 7. L2 norm of the constraints, C , of the polarized 
Gowdy wave tests for the standard BSSN and two adjusted 
formulations. The vertical axis is the logarithm of the L2 
norm of C 2 and the horizontal axis is backward time. The 
solid line (A) is the standard BSSN formulation, the dot- 
ted line (B) is the ^-adjusted BSSN formulation with ka = 
-Id -0 ' 2 , and the dot-dashed line (C) is the C 2 -adjusted BSSN 



formulation with A„ 



-10" 



*a = 



\k = 

14.3 



10" 4ft , A, = -KT 11 , 
W~ 1A , and A F = -10 -1 *" 5 . Note that the signa- 
tures of ka and As are negative since the simulations evolve 
backward. We see that lines (A) and (C) are identical until 
t = —200. Line (C) then decreases and maintains its magni- 
tude under O(10~ 2 ) after t = —400. We confirm this behavior 
until t = -1500. 



by these terms. Therefore, these terms are considered 
to become significant of approximately t = —200 when 
the violations of the A, T~L, and 5-constraints become a 
certain order of magnitude. 

In contrast to the gauge- wave testbed (Figj4j) , we pre- 
pared FiglHl which shows the magnitudes of the ratio of 
the adjusted terms to the original terms. Since the mag- 
nitudes of the adjusted terms of tp and jij can be disre- 
garded, the effect of the reduction of the adjusted terms 
of tp and is negligible. Therefore, the C 2 -adjusted 
BSSN evolution in the Gowdy wave can be regarded as 
maintaining its original hyperbolicity. 

We repeated the parameter-dependency survey of 
(A v , Xk, Xj, A;j, Ap) for this spacetimc evolution. Similar 
to Sec II V A"2l we found that constraint-damping feature 
is sensitive to both Xk and At, of which Xj works ef- 
fectively than Xk- We see the most controlled evolution 



when A^ = 10 1 , 



than that of X^ 



10° or X A = 10" 



FIG. 8. The same with Fig(6]but for the C 2 -adjusted BSSN 
formulation. The parameters, Ax, A^, At, Af), are the 
same with those for (C) in Fig[7] We see that the violation 
of the A4-constraint decreases and becomes the lowest after 



t 



-700. 



-2 

-4 

-6 

-8 
-10 
-12 
-14 



I I I I 
A equation /i A A /\ 


....-^-<p equation ^jjk 






y equation 


r equation 


f. < A, 


1 W * ' 1 





200 



400 600 
-Time 



800 



1000 



FIG. 9. L2 norm of the ratio (adjusted terms)/(original 
terms) of each evolution equation for the C 2 -adjusted BSSN 
formulation, (|3.19I) - (I3.23|) . We see that the largest ratio is 
that for the evolution of Aij. The corrections to the 7^ and 
F 8 evolution equations are reasonably small. 



the violations of all the constraint with cq — ca — cs = 
0. We see that all the violations of the constraints are 
larger than those in FigJSJ This result is consistent with 
the discussion in Scc lIII Bl 



3. Contribution of Algebraic Constraints 
in Definition of C 2 

In Sec lIII Bl we investigated the effect of the defini- 
tion of C 2 . Similar to the gauge- wave tests in the previ- 
ous subsection, we show the effect of constraint damping 
caused by the algebraic constraints. In FigfTUl we plot 



V. SUMMARY AND DISCUSSION 

To obtain an evolution system robust against the vi- 
olation of constraints, we derived a new set of adjusted 
BSSN equations applying the idea proposed by Fiske [2~H 
which we call a "C 2 -adjusted system." That is, we added 
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H-constraint 
M-constraint 
G-constraint 
A-constraint 
S-constraint 



200 
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-Time 



800 



1000 



FIG. 10. Difference with the definition of C with ca = 
ca — cs — 0. The coefficient parameters, \ v , \k, A^, and 
Ap, are all the same as those for (C) in Fig[7] In comparison 
with FigfS] all the violations of the constraints are larger. 



the functional derivatives of the norm of the constraints, 
C 2 , to the evolution equations [ ([3TT!^ - ([3T2"3"]) ]. We per- 
formed numerical tests in the gauge-wave and Gowdy 
wave spacetimes and confirmed that the violations of con- 
straints decrease as expected, and that longer and accu- 
rate simulation than that of the standard BSSN evolution 
is available. 

The construction of the C 2 -adjusted system is straight- 
forward. However, in BSSN, there are two kinetic con- 
straints and three additional algebraic constraints com- 
pared to the ADM system; thus, the definition of C 2 is 
a matter of concern. By analyzing constraint propaga- 
tion equations, we concluded that C 2 should include all 
the constraints. This was also confirmed by numerical 
tests. The importance of such algebraic constraints sug- 
gests the similar treatment when we apply this idea to 
other formulations of the Einstein equation. 

To evaluate the reduction of the violations of the 
constraints, we also compared evolutions with the A- 

22[. We con- 



adjusted BSSN formulation proposed in 
eluded that the C 2 -adjusted BSSN formulation exhibits 
superior constraint damping to both the standard and 
A-adjusted BSSN formulations. In particular, the life- 
times of the simulations of the C 2 -adjustcd BSSN for- 
mulation in the gauge-wave and Gowdy wave testbeds 
are as ten-times and twice as longer than those of the 
standard BSSN formulation, respectively. 

So far, many trials have been reported to improve 
BSSN formulation (e.g. [22|,[3l[). Recently, for example, 
a conformal-traceless Z4 formulation was proposed with 
its test demonstrations [l7|- Among them, Fig.l of [l?} 
can be compared with our FigJ3] [(B-l) and (C-l)] as the 
same gauge-wave test. The violation of H-constraint in 



C 2 -adjusted evolution looks smaller than that of new Z4 
evolution, but regarding the blow-up time of simulations, 
new Z4 system has advantage. 

Fiske reported the applications of the idea of C 2 - 
adjustment to "linearized" ADM and BSSN formulations 
in his dissertation (25|. (As he mentioned, his BSSN is 
not derived from the standard BSSN equations but from 
a linearized ADM using a new variable, V. His set of 
BSSN equations also docs not include the A- and S- 
constraints in our notation.). He observed damping of 
the constraint violation of five orders of magnitude and 
the equivalent solution errors in his numerical evolution 
tests. Our studies show that the full BSSN set of equa- 
tions with fully adjusted terms also produces the desired 
constraint-damping results (Figf2] and FigEJ) , although 
apparent improvements are at fewer orders of ma gnit ude. 

When applied this idea to the ADM system [26|], we 
found that the adjustment to the Kij -evolution equation 
is essential. In the present study, we found that the ad- 
justment to the Ay-evolution equation is essential for 
controlling the constraints. In both cases, the associated 
adjustment parameters (Lagrangian multipliers), A^ in 
this study, are sensitive and require fine-tuning. In fu- 
ture, automatic controlling system such that monitoring 
the order of constraint violations and maintaining them 
by tuning the parameters automatically would be help- 
ful. Applications of control theory in this direction are 
being investigated. 

The correction terms of the C 2 -adjusted system in- 
clude higher-order derivatives and are not quasi-linear; 
thus, little is known mathematically about such systems. 
These additional terms might effectively act as artificial 
viscosity terms in fluid simulations, but might also en- 
hance the violation of errors. To investigate this direction 
further, the next step is to apply the idea to a system in 
which constraints do not include second-order derivatives 
of dynamical variables. We are working on the Kidder- 
Scheel-Teukolsky formulation [To[ as an example of such 
a system, which we will report in the near future. 
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Appendix A: Additional C 2 -adjusted Terms 

Thejadjusted terms_ 5C 2 /Sip, 5C 2 /8K, 5C 2 /5j mn , 
5C 2 /5A rnni and SC 2 /5T a in (I31B - (I3"^31 are written as 
follows: 
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sc 2 

Sip 



sc 2 

~5K 

SC 2 



SC 2 

sc 2 

8T a 



: 2H 1 H - 2{d a H%)% - 2H%d a n + 2(d a d b H% b )H + 2{d a Hf)d b U + 2{d b Hf)d a U + 2Hfd a d b U 

- 2{d a M ll a )e~ iv T J Mj + 8M li a e-^{d aV )^M j - 2M ll a e~^{d a ^)M, - 2M ll a e-^d a M ] 

- l-'-'r ; > M,M, + 4c G e 4v 7 y W, (Al) 

: 2H 4 U - 2{d l M 2l t )e- iip ^M j + 8M 2 /e^(a^)7 !J Mj - 2M- ll l e~^{da ir )M j - 2M 2 /e- 4v j ij d e Mj, 

(A2) 

: 2H™ n H - 2(d i H$ 7in )H - 2Hl mn d t H + 2(d l d 3 H l i > mn )U + 2{d l H l i > mn )d J 'H + 2{d J H l 7 3mn )d t H 

+ 2H % i > mn d l d J U + 2M 3t mn e~ 4v ^M J - 2(d c M il cmn )e- 4 ^ J Mj + 8M At cmn e- Av {d c ip)^ M 3 

- 2M il cmn e- A ' p (d c ^)M 3 - 2M u cmn e- iv Y J dcM 3 - e'^ff M,^ + 2c G G\ mn e 4v % Q 3 

- 2c G {d t GT nl )e Av lvG 1 - 8c G G| w e 4 ^(9^)7ii^ - 2c G G^ nl (da^Q 3 - 2c G G l 2 mnl l e A ^ d& 

+ c G e 4v G m g n + 2c A A™ n A + 2c s S™ n S, (A3) 

: 2H™ n H + 2e~ iip j l3 M 5l mn M J - 2(d c M 6l cmn )e- 4v J 13 Mj + 8M m cmn e- 4v (d c <p)j 13 Mj 

- 2M 6l cmn e- Al(> {d c ^ 3 )M ] - 2M 6l cmn e- 4 ^d c Mj + 2c A A 2 nn A, (A4) 
2H 9a H - 2(d b H b 10a )H - 2H b Wa d b H + 2c G G\ a e^y l3 Q\ (A5) 



where 
Hi 



-Ae- ilp R + 32e~ 4v {D l Diip + (£> i¥ >)(£>V)}, 

(A6) 



iJ 8 mn = -2A mn - (2/3)7 m "K, 
H 9a = (l/2)e- 4 ^7 y - a , 



(A13) 
(A14) 



£T 2 a = 8e- 4v (7 ij r a y - 2D a <p), 



(A7) 



H b = e~^5 h 



(A15) 



i •> 



(A8) 



M u a = 6A a , - 2A mn j mn 5 a . l , (A16) 



Hi = (4/3)# - {2/$)Y 3 A 



i j ■■ 



(A9) 



M 2i j = -(2/3)^, 



(A17) 



- 2e- 4ip f km :j f jn k - 2e~ iv f il{m f n) u 

" 1 1 ai & J- J- li 

+ {l/2)e- A ^ a a l] l am l ln + 8e- 4 ^D m D n V 
+ 2A mb A n b + (2/3) A mn K, 



jjhnn = e -4 V |p^mn + 2 fM< + (l/2)r £ 7 r 



(A10) 



+ 87 £(m ( J D"V) - 47™DV}, (All) 



3j 



-6(D^ m ip)A n h + 2(D l ip)A mn - D^ m A n \ t 



V7 J 



(A18) 



M 4 r n = -^M" 1 ^ + (l/2)7 mn J 4 c l - (l/2)l nm J c i5 

(A19) 

+ (1/2)7™%, (A20) 



M( .cmn = Tjcimgn)^ 



(A21) 



H ijmn = _(xj2) e -^ mn l l \ (A12) 



(A22) 
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G iaU = _^{b~a)i + (l/ 2 )j ab j U , (A23) 

*'r (A24) 

A ab = _£ab^ ^ A25 ^ 

Af = j ab , (A26) 



d t U 

d t M a 
d t G a 

d t A 
d t S 

and those of the A-adjusted BSSN formulation are 



d t H = [Original Terms], (B6) 

d t M t = [Original Terms] + (1/2)k a AM„ (B7) 

d t G i = [Original Terms], (B8) 

d t A = [Original Terms] + n A S ij d l M J , (B9) 

d t S = [Original Terms], (BIO) 

where A is the Laplacian operator in fiat space. "Origi- 



nal Terms" refers to the right-hand side of the constraint 



Sf = (l/2)e a ^ k e bn % n j ki . (A27) 

Appendix B: Constraint Propagation Equations of 
Adjusted BSSN Formulations 

Here we give the constraint propagation equations for 
the C 2 -adjusted BSSN formulation and the A-adjustcd 
BSSN formulation in Minkowskii spacetime. For sim- 
plicity, we Set Xjijmn = Xj5i m 5j n , ^Aijmn = ^A^ im ^3 n ' 

and A~ = Xf,5 tJ . The constraint propagation equations 
of the C 2 -adjustcd BSSN formulation arc 



I 

propagation equations for the standard BSSN formula- 
tion. Full expressions for the terms are given in the ap- 
pendix of |22|. 

Appendix C: Constraint Propagation Equations of 
Standard BSSN Formulation with = 

The constraint propagation equations for the standard 
BSSN formulation with fH 1 = are as follows (the full 
expressions are available in the appendix of [22[). 



[Original Terms] + (-128A V A 2 - (3/2)A^A 2 + 2A F A) H + c G (-(l/2)A^A0 m - 2A F <9 m ) G m + 3c s X^AS, 

(Bl) 

[Original Terms] + ( {8/9)X K S bc d a d b + X x A8 a c + Xj6 bc d a d b \M c - 2c A X x d a A, (B2) 



[Original Terms] + 5 ab ((1/2) X^d b A + 2X f d b ) U + c G (X^A5\ + (l/2)X^S ac d c d b - 2X r S a b ) G b - c s X^S ab d b S, 

(B3) 

[Original Terms] + 2X A S ij {diM 3 ) - 6c A X A A, (B4) 
[Original Terms] + SX^AH + c G X^d t G l - 6c s X^S, (B5) 



dtU = [{2/i)aK + {2/3)aA]n + [-Ae^ 'a(a k ^ kj ~ 2e~ 4 * {d k a)^ k ]Mj 
+ [-2ae-^A k 3 d k - ae~^ '{djA kt )j kt - e-^(d ja )A]G j 

+ fiae-^-^idz^Adk + (l/2) ae - 4 ^7- 1 (a^)7 £fc 5 fc + (\/2)e- i ^- 1 (d l a)^ k Ad k ]S 
+ [(4/9)aKA - {8/9)aK 2 + (4/3)ae~^(d l d J ipW 3 + (8/3)ae- 4v (d k <p){da ek ) + ae'^ (dtf jk )d k 
+ 8ae- A ^ 3k {d^)d k + ae~^ k d 3 d k + ^-^(d^d^)^ + e - iv (d e a)(d k ^ ik ) + 2e~ il f'{d l a)^ k d k 
+ e - 4 ^ ek (dtd k a)]A, (CI) 



12 



d t M t = [-(l/3)(flSia) + (1/6)9^ + aKMi + [ae- A ^ km (d k cp)(d^ mi ) - (l/2)ae-^T m u ^ M (d^ mi ) 
+ (l/2)ae- 4 ^ mk (d k d^ ml ) + (l/2)ae-^- 2 (d l S)(d J S) - {1 / 4)ae~^ {d^ H ){d^ u ) 
+ ae- 4 ^ km (d k <p)^ t d m + ae-^idjtfdi - (l^ae" 4 ^™^^ + ae~^ mk f ijk d m 
+ (\/2)ae- i ^ l % l d k d l + (\/2)e- i ^ mk {d^ lm ){d k a) + {\/2) e -^(d ] a)d l + {\/2)e- i ^ m % l {d k a)d m ]g^ 
+ [-A k i(d k a) + (l/9){ aj )K + {A/9)a(diK) + (l/^aKd, - aA k i d k ]A, (C2) 

d t g l = 2aj ij Mj + [4aj ij (Dj(p) - aFf ij dj - (d k a)^ k ]A, (C3) 
d t A = aKA, (C4) 
d t S = -2ajA. (C5) 
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